Spatial data in R: an introduction to the sf package

---

INSTALLATION AND LOADING

To install sf from CRAN:

install.packages("sf")  

or the development version from GitHub:

devtools::install_github("edzer/sfr")

then load:

library(sf)

NB MacOSX and LINUX users need to install a number of geospatial libraries (GEOS, GDAL, and proj.4).

The tidyverse package also needs to be loaded:

library(tidyverse)

EXAMPLE DATA

A GeoJSON of Greater Manchester’s wards was created from a vector boundary file available from ONS’s Open Geography Portal. The GeoJSON is projected in British National Grid (EPSG:27700) and originally derives from the Ordnance Survey.

Point data are supplied by data.police.uk and represent incidents of anti-social behaviour and crime recorded by the Greater Manchester Police during February 2017. The incidents are supplied with latitude and longitude coordinates which have undergone an anonymisation process.

READING AND WRITING SPATIAL DATA

Reading spatial data (polygons)

bdy <- st_read("data/wards.geojson")

Reading layer `OGRGeoJSON' from data source `/Users/hcpartridge/Desktop/rumgroup/data/wards.geojson' using driver `GeoJSON'
converted into: POLYGON
Simple feature collection with 215 features and 12 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 351664 ymin: 381168.6 xmax: 406087.5 ymax: 421039.8
epsg (SRID):    27700
proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs

bdy

Simple feature collection with 215 features and 12 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 351664 ymin: 381168.6 xmax: 406087.5 ymax: 421039.8
epsg (SRID):    27700
proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs
    objectid    wd16cd                           wd16nm wd16nmw   lad16cd    lad16nm  bng_e  bng_n     long
1        772 E05000851                             Ince    <NA> E08000010      Wigan 359959 405278 -2.60570
2        773 E05000852                       Leigh East    <NA> E08000010      Wigan 367298 400695 -2.49447
3        774 E05000853                      Leigh South    <NA> E08000010      Wigan 366260 398660 -2.50990
4        775 E05000854                       Leigh West    <NA> E08000010      Wigan 363428 400762 -2.55282
5        776 E05000855                      Lowton East    <NA> E08000010      Wigan 362642 397699 -2.56431
6        777 E05000856                           Orrell    <NA> E08000010      Wigan 353321 404083 -2.70568
7        778 E05000857                        Pemberton    <NA> E08000010      Wigan 355101 405469 -2.67903
8        779 E05000858     Shevington with Lower Ground    <NA> E08000010      Wigan 354692 408504 -2.68564
9        780 E05000859           Standish with Langtree    <NA> E08000010      Wigan 355996 410264 -2.66620
10       781 E05000860                        Tyldesley    <NA> E08000010      Wigan 369533 402337 -2.46094
11       782 E05000861                    Wigan Central    <NA> E08000010      Wigan 358550 406484 -2.62712
12       783 E05000862                       Wigan West    <NA> E08000010      Wigan 356532 406901 -2.65763
13       784 E05000863                       Winstanley    <NA> E08000010      Wigan 355407 402755 -2.67402
14       785 E05000864                   Worsley Mesnes    <NA> E08000010      Wigan 357808 403626 -2.63793
15       571 E05000650                    Astley Bridge    <NA> E08000001     Bolton 370670 412906 -2.44479
16       572 E05000651                         Bradshaw    <NA> E08000001     Bolton 374590 413058 -2.38555
17       573 E05000652                       Breightmet    <NA> E08000001     Bolton 374437 409641 -2.38758
18       574 E05000653                    Bromley Cross    <NA> E08000001     Bolton 372113 414150 -2.42309
19       575 E05000654                         Crompton    <NA> E08000001     Bolton 371741 410493 -2.42838
20       576 E05000655                        Farnworth    <NA> E08000001     Bolton 373564 406162 -2.40047
21       577 E05000656                      Great Lever    <NA> E08000001     Bolton 371959 408204 -2.42488
22       578 E05000657                        Halliwell    <NA> E08000001     Bolton 370828 409676 -2.44210
23       579 E05000658                     Harper Green    <NA> E08000001     Bolton 371512 406080 -2.43144
24       580 E05000659               Heaton and Lostock    <NA> E08000001     Bolton 367961 409240 -2.48536
25       581 E05000660             Horwich and Blackrod    <NA> E08000001     Bolton 362140 410545 -2.57343
26       582 E05000661               Horwich North East    <NA> E08000001     Bolton 365750 412038 -2.51906
27       583 E05000662                           Hulton    <NA> E08000001     Bolton 369418 405943 -2.46303
28       584 E05000663                         Kearsley    <NA> E08000001     Bolton 375564 404886 -2.37019
29       585 E05000664     Little Lever and Darcy Lever    <NA> E08000001     Bolton 375274 407414 -2.37476
30       586 E05000665                         Rumworth    <NA> E08000001     Bolton 370267 408025 -2.45041
31       587 E05000666                        Smithills    <NA> E08000001     Bolton 369230 411951 -2.46646
32       588 E05000667            Tonge with the Haulgh    <NA> E08000001     Bolton 372953 409770 -2.41001
33       589 E05000668 Westhoughton North and Chew Moor    <NA> E08000001     Bolton 365620 406875 -2.52045
34       590 E05000669               Westhoughton South    <NA> E08000001     Bolton 364395 404864 -2.53871
35       591 E05000670                           Besses    <NA> E08000002       Bury 381771 405926 -2.27659
36       592 E05000671                           Church    <NA> E08000002       Bury 378564 410535 -2.32531
37       593 E05000672                             East    <NA> E08000002       Bury 382701 411330 -2.26286
38       594 E05000673                            Elton    <NA> E08000002       Bury 379386 412322 -2.31301
39       595 E05000674                         Holyrood    <NA> E08000002       Bury 383274 405503 -2.25388
40       596 E05000675                         Moorside    <NA> E08000002       Bury 381103 412374 -2.28706
41       597 E05000676                      North Manor    <NA> E08000002       Bury 378620 414508 -2.32474
42       598 E05000677                  Pilkington Park    <NA> E08000002       Bury 379384 405697 -2.31260
43       599 E05000678                   Radcliffe East    <NA> E08000002       Bury 378655 408195 -2.32378
44       600 E05000679                  Radcliffe North    <NA> E08000002       Bury 376513 409386 -2.35621
45       601 E05000680                   Radcliffe West    <NA> E08000002       Bury 377876 406692 -2.33543
46       602 E05000681                       Ramsbottom    <NA> E08000002       Bury 378147 417005 -2.33207
47       603 E05000682                         Redvales    <NA> E08000002       Bury 380506 409255 -2.29589
48       604 E05000683                        St Mary's    <NA> E08000002       Bury 380107 403393 -2.30154
49       605 E05000684                          Sedgley    <NA> E08000002       Bury 382536 403044 -2.26488
50       606 E05000685                       Tottington    <NA> E08000002       Bury 377233 412927 -2.34559
51       607 E05000686                         Unsworth    <NA> E08000002       Bury 382436 407933 -2.26667
52       608 E05000687              Ancoats and Clayton    <NA> E08000003 Manchester 387636 398888 -2.18779
53       609 E05000688                          Ardwick    <NA> E08000003 Manchester 385839 396661 -2.21477
54       610 E05000689                          Baguley    <NA> E08000003 Manchester 381171 388024 -2.28457
55       611 E05000690                         Bradford    <NA> E08000003 Manchester 387172 398057 -2.19475
56       612 E05000691                       Brooklands    <NA> E08000003 Manchester 380793 389733 -2.29036
57       613 E05000692                          Burnage    <NA> E08000003 Manchester 386408 392483 -2.20602
58       614 E05000693                      Charlestown    <NA> E08000003 Manchester 387385 403256 -2.19175
59       615 E05000694                         Cheetham    <NA> E08000003 Manchester 384408 400457 -2.23652
60       616 E05000695                         Chorlton    <NA> E08000003 Manchester 381131 393791 -2.28552
61       617 E05000696                    Chorlton Park    <NA> E08000003 Manchester 382585 392536 -2.26356
62       618 E05000697                      City Centre    <NA> E08000003 Manchester 383942 398119 -2.24342
63       619 E05000698                        Crumpsall    <NA> E08000003 Manchester 384444 402511 -2.23608
64       620 E05000699                    Didsbury East    <NA> E08000003 Manchester 384893 390817 -2.22874
65       621 E05000700                    Didsbury West    <NA> E08000003 Manchester 383755 391191 -2.24588
66       622 E05000701                      Fallowfield    <NA> E08000003 Manchester 384495 394300 -2.23490
67       623 E05000702                     Gorton North    <NA> E08000003 Manchester 388555 396421 -2.17385
68       624 E05000703                     Gorton South    <NA> E08000003 Manchester 388320 395173 -2.17734
69       625 E05000704                        Harpurhey    <NA> E08000003 Manchester 386364 401341 -2.20707
70       626 E05000705                  Higher Blackley    <NA> E08000003 Manchester 384361 404086 -2.23741
71       627 E05000706                            Hulme    <NA> E08000003 Manchester 383630 396698 -2.24805
72       628 E05000707                      Levenshulme    <NA> E08000003 Manchester 387196 393796 -2.19421
73       629 E05000708                        Longsight    <NA> E08000003 Manchester 387026 395294 -2.19684
74       630 E05000709  Miles Platting and Newton Heath    <NA> E08000003 Manchester 387660 400026 -2.18748
75       631 E05000710                        Moss Side    <NA> E08000003 Manchester 383884 395430 -2.24416
76       632 E05000711                           Moston    <NA> E08000003 Manchester 387958 401864 -2.18305
         lat st_areasha st_lengths                       geometry
1   53.54262    4780173  11442.625 POLYGON((359198.460435328 4...
2   53.50194    4972092  11056.203 POLYGON((367700.613520717 4...
3   53.48358    6521214  12377.967 POLYGON((364489.980933077 3...
4   53.50228    7315344  16159.924 POLYGON((363107.910792193 3...
5   53.47470   11185415  17990.730 POLYGON((361792.312268886 3...
6   53.53133    8884317  19650.851 POLYGON((353418.76434581 40...
7   53.54394    6239016  12622.229 POLYGON((355533.978524123 4...
8   53.57118    8963755  17516.213 POLYGON((353418.76434581 40...
9   53.58711   10401896  16021.993 POLYGON((356796.803981366 4...
10  53.51683    5559650  12889.007 POLYGON((368822.197554502 4...
11  53.55335    3525250   9567.055 POLYGON((357695.08368452 40...
12  53.55693    2708373   9178.338 POLYGON((356255.095937935 4...
13  53.51957    6183791  10253.866 POLYGON((355533.978524123 4...
14  53.52760    5930122  11368.272 POLYGON((358614.66424552 40...
15  53.61189    6601953  14566.948 POLYGON((370549.775116555 4...
16  53.61346    9283676  16445.880 POLYGON((375026.364570552 4...
17  53.58274    3704010   9695.609 POLYGON((373586.218118696 4...
18  53.62315    6975068  15981.476 POLYGON((373169.622402638 4...
19  53.59026    3483490  12335.029 POLYGON((372158.45537916 40...
20  53.55143    4355269  11771.147 POLYGON((374784.205484859 4...
21  53.56970    4036982  11531.263 POLYGON((370802.588433939 4...
22  53.58287    2820871   7928.388 POLYGON((369793.9063342 409...
23  53.55058    3158073   9727.627 POLYGON((370505.503263172 4...
24  53.57878    8065909  17328.300 POLYGON((366456.699782229 4...
25  53.59012   14745198  24960.144 POLYGON((363194.780215914 4...
26  53.60379   10174487  18821.484 POLYGON((366456.699782229 4...
27  53.54923   10475447  14821.862 POLYGON((368584.593338779 4...
28  53.54006    7692125  14032.845 POLYGON((376604.682156536 4...
29  53.56276    5163985  13467.401 POLYGON((374784.205484859 4...
30  53.56800    1937048   6183.851 POLYGON((370505.503263172 4...
31  53.60322   10977728  20306.973 POLYGON((366341.544716241 4...
32  53.58383    3078336   8853.175 POLYGON((372158.45537916 40...
33  53.55738   16751948  24438.248 POLYGON((363114.143528574 4...
34  53.53922    6313599  13102.652 POLYGON((363400.546355882 4...
35  53.54966    2544591   7588.175 POLYGON((380940.289808304 4...
36  53.59096    4036063  12604.877 POLYGON((375857.068369663 4...
37  53.59826    4195185  13450.955 POLYGON((379584.387973824 4...
38  53.60706    3316545   9011.894 POLYGON((378920.194781534 4...
39  53.54591    4460157  12444.692 POLYGON((382464.27927018 40...
40  53.60759    4703167  10772.104 POLYGON((381389.068949916 4...
41  53.62667    9789443  23196.083 POLYGON((375026.364570552 4...
42  53.54751    7217589  15340.667 POLYGON((380774.396924753 4...
43  53.56993    6470312  14319.963 POLYGON((379233.131869085 4...
44  53.58055    7053574  13177.442 POLYGON((375332.182010675 4...
45  53.55639    4896161  12607.624 POLYGON((377289.29754458 40...
46  53.64910   13489427  21951.402 POLYGON((375569.530811583 4...
47  53.57953    3676876  11852.667 POLYGON((381691.379115022 4...
48  53.52683    5014610  11047.551 POLYGON((380940.289808304 4...
49  53.52378    2110743   8605.973 POLYGON((381671.176295511 4...
50  53.61241    7792259  17558.920 POLYGON((375857.068369663 4...
51  53.56772    8696448  16181.589 POLYGON((381465.482109123 4...
52  53.48657    4434567  15297.173 POLYGON((388979.59511252 39...
53  53.46651    4090129   8290.054 POLYGON((385048.700403615 3...
54  53.38873    3819630   9646.251 POLYGON((382222.444481428 3...
55  53.47909    5160234  12321.011 POLYGON((389082.893464996 3...
56  53.40408    4166786  13208.062 POLYGON((379868.483208293 3...
57  53.42897    2324455   7870.979 POLYGON((385852.697738632 3...
58  53.52583    3572634   8915.117 POLYGON((388814.218263555 4...
59  53.50059    4703758  10305.930 POLYGON((385383.442175553 3...
60  53.44056    2495801   8555.880 POLYGON((382069.783909148 3...
61  53.42933    4798682  11399.285 POLYGON((382951.464352811 3...
62  53.47956    3102193   8569.712 POLYGON((384111.586234225 3...
63  53.51905    3186206  10323.461 POLYGON((383059.682257273 4...
64  53.41395    3542769  10341.897 POLYGON((384160.163612966 3...
65  53.41728    2823282  10970.634 POLYGON((384132.27552782 39...
66  53.44525    2265053   8571.091 POLYGON((382951.464352811 3...
67  53.46442    3520488   9673.390 POLYGON((390196.360627751 3...
68  53.45319    3315488  11618.528 POLYGON((387559.543108503 3...
69  53.50859    3815368   9202.353 POLYGON((385383.442175553 3...
70  53.53320    7040635  13779.614 POLYGON((382464.27927018 40...
71  53.46678    2396423   7806.111 POLYGON((383254.103669637 3...
72  53.44079    2356776   9676.619 POLYGON((386989.75973683 39...
73  53.45425    1498778   6469.738 POLYGON((385590.186839664 3...
74  53.49680    4353329  12542.810 POLYGON((384894.453135228 3...
75  53.45539    1825331   7563.015 POLYGON((383254.103669637 3...
76  53.51333    3093419  10756.913 POLYGON((388794.32215659 40...
 [ reached getOption("max.print") -- omitted 139 rows ]

Reading data with coordinates (points)

pts <- read_csv("data/2017-03-greater-manchester-street.csv")
pts <- st_as_sf(pts, coords = c("Longitude", "Latitude"))
pts

Simple feature collection with 35094 features and 10 fields
geometry type:  POINT
dimension:      XY
bbox:           xmin: -2.716562 ymin: 53.34144 xmax: -1.96473 ymax: 53.66877
epsg (SRID):    NA
proj4string:    NA

Writing spatial data

st_write(bdy, "boundaries.shp")

SF OBJECTS

Attribute table (dataframe) AND geometry (list-column) with coordinates, CRS, bbox

st_geometry(bdy) # print geometry

Geometry set for 215 features 
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 351664 ymin: 381168.6 xmax: 406087.5 ymax: 421039.8
epsg (SRID):    27700
proj4string:    +proj=tmerc +lat_0=49 +lon_0=-2 +k=0.9996012717 +x_0=400000 +y_0=-100000 +ellps=airy +towgs84=446.448,-125.157,542.06,0.15,0.247,0.842,-20.489 +units=m +no_defs
First 5 geometries:

MANIPULATE USING DPLYR

Rename variables

glimpse(bdy)

Observations: 215
Variables: 4
$ ward        <fctr> Ince, Leigh East, Leigh South, Leigh West, Lowton East, Orrell, Pemberton, Shevington wi...
$ census_code <fctr> E05000851, E05000852, E05000853, E05000854, E05000855, E05000856, E05000857, E05000858, ...
$ borough     <fctr> Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wigan, Wiga...
$ geometry    <simple_feature> POLYGON((359198.460435328 4..., POLYGON((367700.613520717 4..., POLYGON((36448...

Count frequency of wards by borough

bdy %>% 
  group_by(borough) %>% 
  count() %>% 
  arrange(desc(n)) %>% # sort in descending order
  st_set_geometry(., NULL) # hide geometry

Select features

chorlton <- bdy %>% 
  filter(ward == "Chorlton")
plot(st_geometry(bdy))
plot(st_geometry(chorlton), col = "red", add = TRUE)

Using sf functions to add geometry column to dplyr chain

bdy <- bdy %>% 
  mutate(area = st_area(.)) # returns the area of a feature
bdy %>% 
  select(ward, area) %>% 
  arrange(desc(area)) %>% 
  slice(1:10) %>% 
  st_set_geometry(., NULL)

PROJECTION

Check and assign CRS

st_crs(pts) 

$epsg
[1] NA

$proj4string
[1] NA

attr(,"class")
[1] "crs"

pts <- st_set_crs(pts, 4326) # assign Lat/Long (epsg:4326)
st_crs(pts)

$epsg
[1] 4326

$proj4string
[1] "+proj=longlat +datum=WGS84 +no_defs"

attr(,"class")
[1] "crs"

Reproject CRS

bdy_WGS84 <- st_transform(bdy, 4326)
st_crs(bdy_WGS84)

$epsg
[1] 4326

$proj4string
[1] "+proj=longlat +datum=WGS84 +no_defs"

attr(,"class")
[1] "crs"

CONVERT TO AND FROM SP OBJECTS

Convert to and from sp objects

bdy_sp <- as(bdy, 'Spatial')
class(bdy_sp)

[1] "SpatialPolygonsDataFrame"
attr(,"package")
[1] "sp"

bdy_sf <- st_as_sf(bdy_sp)
class(bdy_sf)

SPATIAL OPERATIONS

Buffer features

buffer <- chorlton %>% 
  st_buffer(dist = 1000)
plot(st_geometry(buffer))
plot(st_geometry(chorlton), col = "red", add = TRUE)

Buffer and intersect

pts_sub <- bdy %>%
  filter(ward == "Chorlton") %>%
  st_buffer(dist = 1000) %>%
  st_intersection(st_transform(pts, 27700)) # reproject pts to BNG (epsg:27700)
plot(st_geometry(buffer))
plot(st_geometry(chorlton), col = "red", add = TRUE)

plot(st_geometry(pts_sub), col = "black", add = TRUE)

pts_sub %>% 
  group_by(`Crime type`) %>%
  count() %>% 
  arrange(desc(n)) %>% 
  st_set_geometry(., NULL)

Points in polygon

pts %>% 
  filter(`Crime type` == "Vehicle crime") %>%  
  st_join(bdy_WGS84, ., left = FALSE) %>% 
  count(ward) %>% 
  arrange(desc(n)) %>% 
  st_set_geometry(., NULL)

bdy_pts <- pts %>% 
  filter(`Crime type` == "Vehicle crime") %>%  
  st_join(bdy_WGS84, ., left = FALSE) %>% 
  count(ward)

PLOTTING

Base plots

plot(bdy_pts) # plots small multiples if dataframe has several attributes

plot(bdy_pts["n"]) # select the appropriate attribute to plot a single map

library(RColorBrewer) ; library(classInt)
pal <- brewer.pal(5, "RdPu")
classes <- classIntervals(bdy_pts$n, n=5, style="pretty")$brks
plot(bdy_pts["n"], 
     col = pal[findInterval(bdy_pts$n, classes, all.inside=TRUE)], 
     main = "Vehicle crime in Greater Manchester\nMarch 2017", axes = F)
legend("bottomright", legend = paste("<", round(classes[-1])), fill = pal, cex = 0.7) 

ggplot2

devtools::install_github("tidyverse/ggplot2") # NB need development version for geom_sf()

ggplot(bdy_pts) +
  geom_sf(aes(fill = n)) +
  scale_fill_gradientn('Frequency', colours=RColorBrewer::brewer.pal(5,"RdPu"), 
                       breaks = scales::pretty_breaks(n = 5)) +
  labs(fill = "Frequency",
       title = "Vehicle crime",
       subtitle = "March 2017",
       caption = "Contains OS data © Crown copyright and database right (2017)") +
  theme_void()

plotly

library(plotly)
p <- ggplot(bdy_pts) +
  geom_sf(aes(fill = n, text = paste0(ward, "\n", "Crimes: ", n))) +
  scale_fill_gradientn('Frequency', colours=RColorBrewer::brewer.pal(5,"RdPu"), 
                       breaks = scales::pretty_breaks(n = 5)) +
  labs(fill = "Frequency",
       title = "Vehicle crime",
       subtitle = "March 2017",
       caption = "Contains OS data © Crown copyright and database right (2017)") +
  theme_void() + 
  coord_fixed(1.3)
ggplotly(p, tooltip = "text")

leaflet

library(leaflet)
pal <- colorBin("RdPu", domain = bdy_pts$n, bins = 5, pretty = TRUE)
leaflet(bdy_pts) %>% 
  addTiles(urlTemplate = "http://{s}.tiles.wmflabs.org/bw-mapnik/{z}/{x}/{y}.png",
    attribution = '&copy; <a href="http://www.openstreetmap.org/copyright">OpenStreetMap</a>, <a href="https://www.ons.gov.uk/methodology/geography/licences">Contains OS data © Crown copyright and database right (2017)</a>') %>% 
  addPolygons(fillColor = ~pal(n), fillOpacity = 0.8,
              weight = 1, opacity = 1, color = "black",
              label = ~as.character(ward)) %>% 
  addLegend(pal = pal, values = ~n, opacity = 0.7, 
            title = 'Vehicle crime (March 2017)', position = "bottomleft") %>%
  addMiniMap(tiles = providers$CartoDB.Positron, toggleDisplay = TRUE)